The birhythmicity increases the diversity of p53 oscillation induced by DNA damage
Wang Dao-Guang1, 2, Zhou Chun-Hong2, Zhang Xiao-Peng1, 3, †
National Laboratory of Solid State Microstructure and Department of Physics, Nanjing University, Nanjing 210093, China
School of Physics and Electronic Engineering, Jiangsu Normal University, Xuzhou 221116, China
Kuang Yaming Honors School, Nanjing University, Nanjing 210023, China

 

† Corresponding author. E-mail: zhangxp@nju.edu.cn

Abstract

The tumor suppressor p53 mediates the cellular response to various stresses. It was experimentally shown that the concentration of p53 can show oscillations with short or long periods upon DNA damage. The underlying mechanism for this phenomenon is still not fully understood. Here, we construct a network model comprising the ATM-p53-Wip1 and p53-Mdm2 negative feedback loops and ATM autoactivation. We recapitulate the typical features of p53 oscillations including p53 birhythmicity. We show the dependence of p53 birhythmicity on various factors such as the phosphorylation status of ATM. We also perform stochastic simulation and find the noise-induced transitions between two modes of p53 oscillation, which increases the p53 variability in both the amplitude and period. These results suggest that p53 birhythmicity enhances the responsiveness of p53 network, which may facilitate its tumor suppressive function.

1. Introduction

Periodic oscillations are involved in various cellular processes such as glucose metabolism,[1] cell cycle,[2] and circadian/circannual rhythm,[3] with the period ranging from less than a second to years. Typically, oscillations have a relatively fixed period under the same external condition. Nevertheless, there may coexist two different stable oscillations, which alternate temporally or one of which predominates depending on the initial condition. This phenomenon is called birhythmicity,[4,5] resulting from the association of bistability and oscillatory dynamics. Birhythmicity has been observed in a number of biochemical systems[6] as well as the heart and nervous systems.[7,8]

The p53 protein is one of the most important tumor suppressors. It can mediate the cellular response to a wide variety of stress signals including DNA damage and oncogene activation.[911] It is increasingly evident that the dynamics of p53 levels play a key role in cell-fate decision. In unstressed cells, p53 is kept at low levels because of its negative regulator Mdm2; p53 induces the production of Mdm2, whereas Mdm2 targets p53 for degradation. Upon DNA damage, the concentrations of p53 and Mdm2 can exhibit a series of discrete pulses in diverse cell lines such as MCF-7, A549 and H1299 ones.[4,12] p53 oscillations also exhibit variability in both the amplitude and period.[12,13] Strikingly, two oscillatory modes have been observed simultaneously: a low-frequency pulsing with the period of ∼ 10 h and a high-frequency pulsing with the period of ∼ 6 h.[13] The underlying mechanism and functional implications of such p53 birhythmicity have attracted much attention.[1417]

Theoretical models were built to interpret the experimental obervations.[14,15] In Ref. [14], the coupling between the positive and negative feedback loops engaging p53 and Mdm2 were taken into account. But it is well known that the ataxia mutant (ATM) kinase plays a critical role in the DNA damage response. Upon DNA damage, inactive ATM dimer converts to active, phosphorylated monomer via intermolecular auto-phosphorylation (i.e., engaging a positive feedback loop (PFL));[18] ATM then phosphorylates p53 and Mdm2, enhancing the stability and transcriptional activity of p53. Activated p53 induces the production of the Wip1 phosphatase, which in turn dephosphorylates p53 and ATM. It was confirmed that the ATM-p53-Wip1 negative feedback loop (NFL) is indispensable for p53 oscillations.[19,20] Thus, it is essential to explore p53 oscillation within a wider context including ATM and Wip1.

Here, we build a model of p53 network in response to DNA damage, comprising the ATM-p53-Wip1 and p53-Mdm2 NFLs and ATM autoactivation modules. We probe the underlying mechanism for birhythmicity in p53 dynamics and identify the influence of ATM and Wip1. Stochastic simulations show that the birhythmicity in p53 oscillations are helpful for the reproduction of remarkable p53 variability. The birhythmicity makes it possible for p53 levels to switch quickly from one mode to the other following different DNA damage, enhancing the flexibility of the p53 response. This is of functional significance for p53 to modulate the cellular response.

2. Model and method
2.1. Model

Based on previous studies,[13,2022] we construct a model comprising core components associated with p53 oscillation [Fig. 1]. In unstressed cells, p53 remains at low levels due to the negative regulation by the E3 ubiquitin ligase Mdm2.[23] Upon ionizing radiation (IR) or DNA-damaging agents like neocarzinostatin (NCS),[12] DNA double-strand breaks (DSBs) are generated. DNA repair proteins are recruited to fix the damage. Meanwhile, rapid autophosphorylation of ATM is induced,[18,24] resulting in the conversion from ATM dimer (ATM2) to active phosphorylated ATM monomer (ATMp).[24,25] ATMp phosphorylates p53 and Mdm2,[25] accelerating the degradation of nuclear Mdm2 and enhancing the stability and transcriptional activity of p53.[26] Phosphorylated p53 (p53p) induces the production of Mdm2 and Wip1,[27] and Wip1 dephosphorylates p53 and ATMp. Taken together, there are two NFLs and one PFL in the model. Since it was shown experimentally that p53 oscillation can persist for a long time,[13] we suspect that the repair machinery therein may be defective, and thus we do not consider the repair process and take the number of initial DSBs as the input.

Fig. 1. (color online) Schematic diagram of the network model. DNA damage activates ATM, which phosphorylates and accelerates degradation of nuclear Mdm2, thus stabilizing and activating p53. Phosphorylated p53 induces the production of Mdm2 and Wip1, and Wip1 dephosphorylates ATM and p53.
2.2. Deterministic equations

The dynamics of the model are governed by the following ordinary differential equations (ODEs). The kinetic equations for ATM activation are as follows: where [ ] denotes the concentration of components and the function is defined as . Since the total amount of ATM does not change much after IR,[28,29] we assume [ATM]tot = [ATM]+[ATMp]+2[ATM2] as a constant. The binding and disassociation reactions are described according to mass action law, while the enzymatic reactions are characterized by Michaelie–Menten (MM) kinetics. Upon DNA damage, ATM dimers are recruited to damaged DNA and dissociate into monomers,[25,30] which is characterized by the first term on the right hand side of Eq. (1), with denoting the number of DSBs.

The dynamics of p53 are governed by These equations describe the basal transcription, translation, degradation and reversible phosphorylation of p53.

For Mdm2, the equations read The function is defined as . We consider three forms of Mdm2: cytoplasmic Mdm2 (Mdm2c), phosphorylated Mdm2c (Mdm2cp) and nuclear Mdm2 (Mdm2 ). The transcription of mdm2 by p53 is characterized by Hill function with n denoting the degree of cooperativity, while the nuclear import and export of Mdm2 are described by linear terms with rate constants and , respectively.

Similarly, the dynamics of Wip1 are governed by The transcription of wip1 by p53 is modeled by Hill equation with n = 3. All the parameter values are listed in Table 1. The time is in units of minutes, while the concentrations are normalized to be dimensionless.

Table 1.

Default values of parameters.

.
3. Results
3.1. Overview of p53 dynamics

Without DSB, i.e., , the system stays in a stable steady state with kept at a very low level. When , undergoes Hopf bifurcation (HB) [Fig. 2(a)]. If , the steady state becomes unstable, and sustained limit-cycle oscillation appears. Both the amplitude and period of p53 oscillation rise with increasing for . Remarkably, for there exist two branches of periodic solutions: one with larger period and amplitude, and the other with smaller period and amplitude. That is, exhibits birhythmicity in this range. For , only undergoes large-amplitude oscillations, and the period drops slowly with increasing .

Fig. 2. (color online) Deterministic analysis of the model. (a) Bifurcation diagram of [ ] versus . The black solid and dashed curves denote the stable and unstable steady states, respectively. The filled circle denotes the Hopf bifurcation point. The blue curves denote the maximum and minimum of [ ] in limit cycles. The inset shows the oscillation period. (b) Time courses of protein concentrations at .

Figure 2(b) shows the time courses of protein concentrations at , where the oscillation period is 324 min (∼ 5.4 h), agreeing with the typical experimental result.[13] ATM is first activated by DNA damage. ATMp promotes the degradation of nuclear Mdm2 and accumulation of p53p. p53 induces Mdm2 and Wip1, and Wip1 dephosphorylates ATM and p53, leading to a drop in [p53p] and [ATMp]. ATM and p53 can be reactivated by residual DNA damage, and a new round of oscillations is triggered. Collectively, the ATM-p53-Wip1 NFL plays a critical role in p53 pulsing.

We further take a closer look at p53 birhythmicity. For , two different oscillatory modes can be induced at the same . The oscillation with smaller amplitude and period is called limit cycle A (LCA), and the oscillation with larger amplitude and period is called limit cycle B (LCB) [Fig. 3(a)]. The ratio of two periods is 0.765, which is around the experimental value (see Fig. S3 in Ref. [13]). In the phase plane of [p53p] versus [ATMp], two isolated closed trajectories[31] can be seen [Fig. 3(b)]. These two stable limit cycles are separated by an unstable limit cycle with different amplitudes and periods. Clearly, the selection of either mode depends on the initial condition of the system state (see Table 2).

Fig. 3. (color online) With different initial values, the p53 network shows distinct oscillatory modes. (a) Temporal evolution of [ ], [ATMp], and [Mdm2n] under different initial conditions. (b) Trajectories in the [p53p]–[ATMp] planes. Two separate limit cycles are shown by the lines in different colors.
Table 2.

Initial values for different limit cycles.

.
3.2. Phosphorylation of ATM greatly affects the birhythmicity

Although it was previously reported that a single NFL (such as the p53-Mdm2 or p53-Wip1 NFL) can allow for birhythmicity,[13,15,16] more robust birhythmicity can be realized in the system of interlinked NFL and PFL. In our model, the birhythmicity appears via the saddle-node bifurcation on the limit cycle, where a single oscillatory domain can be split into two parts. The similar mechanism was reported in the autocatalytic enzymatic reaction model,[4] the OAK model,[14] and the p53-Mdm2 model.[17] In the following, we explore the dependence of birhythmicity on the features of feedback loops, especially on ATM activity.

We first probe the influence of phosphorylation strength of ATM, which is affected by and . As seen in the bifurcation diagram, [p53p] remains at a low level for or at a high level for [Fig. 4(a)]. At intermediate values, [p53p] undergoes periodic oscillation. The amplitude sharply rises to a high level and then changes slightly over a wide range of . Compared with the amplitude, the period depends more heavily on . Specifically, the birhythmicity appears for .

Fig. 4. (color online) Dependence of [p53p] on and with . (a) Bifurcation diagram of [p53p] versus characterizing the strength of ATM phosphorylation. (b) Bifurcation diagram of [p53p] versus representing the strength of ATM dephosphorylation. The inset shows the oscillation period. The symbols ‘HB’ and ‘SN’ denote Hopf bifurcation and saddle-node bifurcation, respectively.

The dephosphorylation rate of ATMp, , also affects the existence of birhythmicity. Similarly, [p53p] is kept in a steady state when takes lower or higher values, and shows oscillations via the saddle-node (SN) bifurcation. The amplitude and period first drop and then change slightly with increasing until 0.729. The birhythmicity exists for . Together, regulating the phosphorylation strength of ATM can influence the existence of birhythmicity.

We also show the dependence of oscillation period on for different . For low phosphorylation rate (e.g., ), the shorter-period oscillation regime L1, the birhythmicity regime L2, and the longer-period region L3 appear successively with increasing [Fig. 5(a)]. For , there first appears the L3 regime, and the birhythmicity regime moves rightward and widens [Fig. 5(b)]. For , the birhythmicity regime narrows, and the L1 regime occurs at relatively high values [Fig. 5(c)]. For high phosphorylation rate ( ), the birhythmicity regime disappears and there exists only the L1 regime [Fig. 5(d)]. Collectively, the birhythmicity regime L2 is sensitive to the phosphorylation status of ATM.

Fig. 5. (color online) Dependence of the oscillation period on for different at . (a) , (b) , (c) , (d) .

It is worth noting that also affects the presence of birhythmicity. As seen in the two-dimensional phase diagrams in Fig. 6, [p53p] may remain in a steady state (I), an excitable state (II), or oscillatory state (III), depending on the combination of and . p53 oscillations can exist in a relatively wide regime. In the diagram of versus , [p53p] can also remain in a bistable state besides the above three states. On the other hand, in the diagram of versus (with ), the oscillation state can exist in a wide regime, including the birhythmicity regime [Fig. 6(c)]. For , a similar diagram is obtained [Fig. 6(d)]. The birhythmicity regime is close to the boundary between the oscillation and steady states, which suggests that the birhythmicity originates from the SN bifurcation on the limit cycle. This is different from the conclusion in previous studies,[4,6] where one NFL and one PFL are coupled together and the birhythmicity regime is close to the excitability; such a difference is related to distinct network architecture. In our model, the ATM-p53-Wip1 and p53-Mdm2 NFLs are coupled with the PFL involving ATM activation, such that both robust oscillations and birhythmicity can be induced.

Fig. 6. (color online) Analysis of the system state in two-dimensional bifurcation diagrams. (a) - , (b) , (c)–(d) with and 40. The symbols ‘I, II, III, IV’ denote the parameter region corresponding to monostability, excitability, oscillation, and bistability respectively. ‘III-a’ denotes the parameter region corresponding to oscillation with birhythmicity.
3.3. Noise-induced transition between two limit cycles

We have probed the network dynamics in the deterministic case. Given noise always exists in biological systems, it is essential to explore the effect of noise on system dynamics. The noise results from low numbers of reactant molecules or/and perturbations such as changes in parameter values. To this end, we rewrite the dynamic equations as follows: where fi is the same as the term on the right-hand side of Eq. (i), and denotes the noise.[3234] It was shown that noises are often colored, i.e., they have a non-zero correlation time, which can remarkably affect the stochasticity in the system dynamics.[33] For simplicity, the noise in [p53p] is assumed to be colored, whereas the noise in the others is Gaussian white. The colored noise obeys the Ornstein-Uhlenbeck (OU) process:[35,36] where is Gaussian white noise, τ is the correlation time,[37] and D is the noise intensity. The Gaussian white noise has zero mean, i.e., and . Each noise is independent of each other (but with the same D). Numerical simulation of these stochastic equations is performed using the method in Ref. [38].

Without loss of generality, we set D = 0.01 and τ = 10 min in simulation. For different , we plot the time evolution of [p53p] and show the histograms for the amplitude and period of oscillations in Fig. 7. Notably, [p53p] can switch between low- and high-amplitude oscillations. At or 50, most oscillations have low and high amplitudes, respectively, whereas at , the oscillations exhibit larger variability, and the large-amplitude oscillations are nearly comparable to low-amplitude ones. For comparison, the distributions of the period change less remarkably. Our results are consistent with the data in Ref. [13]. Together, p53 may exhibit large variability in oscillations in the presence of noise, and the birhymicity may contribute markedly to this variability. This represents a new source of p53 variability.

Fig. 7. Noise-induced variability in p53 oscillations. The first column shows the temporal evolution of [p53p] at (from top to bottom). The second and third columns show the histograms of the amplitude and period, respectively. In each case, 500 runs of simulations lasting 3000 min are performed.
4. Conclusion

In this article, we have investigated the p53 dynamics in response to sustained DNA damage. Because of the ATM-p53-Wip1 NFL and ATM autoactivation, p53 can exhibit complex dynamics, including simple oscillation and birhythmicity. We probed the dependence of p53 birhythmicity on network parameters and found that the birhythmicity depends on the phosphorylation status of ATM. Here the birhythmicity results from the saddle-node bifurcation on the limit cycle,[39] which is different from the mechanism reported in Refs. [14] and [17], where the system is composed of one PFL and one NFL, both involving p53 and Mdm2. Clearly, our model is built on the experimental finding that the ATM-p53-Wip1 NFL plays an essential role in p53 oscillation and p53 oscillation is driven by ATM oscillation.[19] Our results shed new light on p53 birhythmicity.

We also performed stochastic simulation and exhibited p53 variability in both the amplitude and period. Simulation results agree with the experimental observations. p53 birhythmicity increases the variability and enriches the diversity of p53 oscillation, which may have functional significance since p53 dynamics are closely associated with cell-fate decision. It would be interesting to extend the current model by incorporating downstream effectors of p53 and to explore wether and how p53 birhythmicity contributes to its tumor-suppressing function.

Acknowledgment

The authors thank Prof. Feng Liu for constructive suggestions and Dr. Bo Huang for valuable discussions.

Reference
[1] Schisano B Tripathi G McGee K McTernan P G Ceriello A 2011 Diabetologia 54 1219
[2] Kamb A Gruis N A Weaver-Feldhaus J Liu Q Y Harshman K Tavtigian S V Stockert E Day R S Johnson B E Skolnick M H 1994 Science 264 436
[3] Lincoln G A Clarke I J Hut R A Hazlerigg D G 2006 Science 314 1941
[4] Moran F Goldbeter A 1984 Biophys. Chem. 20 149
[5] Leloup J C Goldbeter A 1999 J. Theor. Biol. 198 445
[6] Decroly O Goldbeter A 1982 Proc. Natl. Acad. Sci. USA 79 6917
[7] Haberichter T Marhl M Heinrich R 2001 Biophys. Chem. 90 17
[8] Goldbeter A Morán F 1988 Eur. Biophys. J. 15 277
[9] Kuerbitz S J Plunkett B S Walsh W V Kastan M B 1992 Proc. Natl. Acad. Sci. USA 89 7491
[10] Levine A J 1997 Cell 88 323
[11] Serrano M Lin A W McCurrach M E Beach D Lowe S W 1997 Cell 88 593
[12] Loewer A Karanam K Mock C Lahav G 2013 BMC Biol. 11 114
[13] Geva-Zatorsky N Rosenfeld N Itzkovitz S Milo R Sigal A Dekel E Yarnitzky T Liron Y Polak P Lahav G Alon U 2006 Mol. Syst. Biol. 2 0033
[14] Ouattara D A Abou-Jaoudé W Kaufman M 2010 J. Theor. Biol. 264 1177
[15] Kim J K Jackson T L 2013 PLoS One 8 e65242
[16] Geva-Zatorsky N Dekel E Batchelor E Lahav G Alon U 2010 Proc. Natl. Acad. Sci. USA 107 13550
[17] Abou-Jaoudé W Chaves M Gouzé J L 2011 PLoS One 6 e17075
[18] Bakkenist C J Kastan M B 2003 Nature 421 499
[19] Batchelor E Mock C S Bhan I Loewer A Lahav G 2008 Mol. Cell 30 277
[20] Zhang X P Liu F Wang W 2011 Proc. Natl. Acad. Sci. USA 108 8990
[21] Batchelor E Loewer A Mock C Lahav G 2011 Mol. Syst. Biol. 7 488
[22] Zhang X P Liu F Cheng Z Wang W 2009 Proc. Natl. Acad. Sci. USA 106 12245
[23] Vassilev L T Vu B T Graves B Carvajal D Podlaski F Filipovic Z Kong N Kammlott U Lukacs C Klein C Fotouhi N Liu E A 2004 Science 303 844
[24] Kozlov S V Graham M E Jakob B Tobias F Kijas A W Tanuji M Chen P Robinson P J Taucher-Scholz G Suzuki K So S Chen D Lavin M F 2011 J. Biol. Chem. 286 9107
[25] Shiloh Y Ziv Y 2013 Nat. Rev. Mol. Cell Biol. 14 197
[26] Stommel J M Wahl G M 2004 EMBO J. 23 1547
[27] Shreeram S Demidov O N Hee W K Yamaguchi H Onishi N Kek C Timofeev O N Dudgeon C Fornace A J Anderson C W Minami Y Appella E Bulavin D V 2006 Mol. Cell 23 757
[28] Banin S Moyal L Shieh S-Y Taya Y Anderson C W Chessa L Smorodinsky N I Prives C Reiss Y Shiloh Y Ziv Y 1998 Science 281 1674
[29] Canman C E Lim D-S Cimprich K A Taya Y Tamai K Sakaguchi K Appella E Kastan M B Siliciano J D 1998 Science 281 1677
[30] Dupré A Boyer-Chatenet L Gautier J 2006 Nat. Struct. Mol. Biol. 13 451
[31] Wang L F Qiu K Jia Y 2017 Chin. Phys. 26 030503
[32] Elowitz M B Levine A J Siggia E D Swain P S 2002 Science 297 1183
[33] Rosenfeld N Young J W Alon U Swain P S Elowitz M B 2005 Science 307 1962
[34] Bowsher C G Voliotis M Swain P S 2013 PLoS. Comput. Biol. 9 e1002965
[35] Uhlenbeck G E Ornstein L S 1930 Phys. Rev. 36 823
[36] Gillespie D T 1996 Phys. Rev. 54 2084
[37] Li W Zhang M T Zhao J F 2017 Chin. Phys. 26 090501
[38] Fox R F 1991 Phys. Rev. 43 2649
[39] Biswas D Banerjee T Kurths J 2017 Chaos 27 063110